home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 5 / Amiga Tools 5.iso / grafik / converter / limbo4.0 / src / 3d / classify.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-07-16  |  5.4 KB  |  187 lines
  1. #include "includes.h"
  2.  
  3.  unsigned long *ML;
  4.  long double *M;
  5.  long double *CM;
  6.  long double *NCM;
  7.  long double *n;
  8.  
  9.  long double Pow2(long double d)
  10.   {
  11.    return (long double)(d*d);
  12.   }
  13.  
  14. /**************************************|****************************************
  15. Routine   : Classify
  16. Input  par: BitMap *SubImage (pointer to structure which determines image to classify)
  17. Output par: Int (Integer which determines the class of 'SubImage')
  18. Function  : Determines the class of a given subimage. Possible classes: EDGE, SHADE.
  19. ***************************************|***************************************/
  20.  
  21.  float *Classify3D(BitMap3D *Image,int XStart,int YStart,int ZStart, int SubSize,
  22.                    int Dim, float TSE)
  23.   {
  24.    register unsigned int x,y,z;
  25.    unsigned long f,ff=0,xx,yy,xy,xxy,xyy,m=0;
  26.    long double mx,my,iM2,iM3,temp;
  27.    static int FirstCall=TRUE;
  28.    float *Features=GimmeAFloatArray(Dim+PREDEFFEATURES);
  29.    double Volume=SubSize*SubSize*SubSize;
  30.    
  31.    if(FirstCall)
  32.     {
  33.      FirstCall=FALSE;
  34.      ML =GimmeAULongArray(10);
  35.      M  =GimmeADoubleArray(10);
  36.      CM =GimmeADoubleArray(10);
  37.      NCM=GimmeADoubleArray(7);
  38.      n  =GimmeADoubleArray(7);
  39.     }
  40.    
  41.    for(x=0;x<10;x++) ML[x]=0;
  42.    
  43.    for (x=0;x<SubSize;x++)
  44.    for (y=0;y<SubSize;y++)
  45.     {
  46.      ML[0] += Image->Map[x+XStart][y+YStart][ZStart]; /* m_00 */      
  47.      for (z=0;z<SubSize;z++) 
  48.       {
  49.        f=Image->Map[x+XStart][y+YStart][z+ZStart];
  50.        ff += f*f;
  51.        m  += f;
  52.       }
  53.     }
  54.    
  55.    Features[MEAN]=(float)m/Volume;
  56.    Features[VAR] =((float)ff-Volume*Features[MEAN]*Features[MEAN])/(Volume-1.0);
  57.    
  58.    if (Features[VAR]>TSE) /* is this an edge ? */
  59.     {
  60.      for (x=0;x<SubSize;x++)
  61.      for (y=0;y<SubSize;y++)
  62.      /* for (z=0;z<SubSize;z++) */
  63.       {
  64.        f=Image->Map[x+XStart][y+YStart][ZStart];
  65.        xx=(x+1)*(x+1);
  66.        yy=(y+1)*(y+1);
  67.        xy=(x+1)*(y+1);
  68.        xxy=xx*(y+1);
  69.        xyy=(x+1)*yy;
  70.        ML[0] += f;            /* m_00 */
  71.        ML[1] += (x+1)*f;      /* m_10 */
  72.        ML[2] += (y+1)*f;      /* m_01 */
  73.        ML[3] += (x+1)*(y+1)*f;/* m_11 */
  74.        ML[4] += xx*f;         /* m_20 */
  75.        ML[5] += yy*f;         /* m_02 */
  76.        ML[6] += xxy*f;        /* m_21 */
  77.        ML[7] += xyy*f;        /* m_12 */
  78.        ML[8] += (x+1)*xx*f;   /* m_30 */
  79.        ML[9] += (y+1)*yy*f;   /* m_03 */
  80.       }
  81.      
  82.      Features[STDDEV]=sqrt(Features[VAR]);
  83.      
  84.      for(x=0;x<10;x++) M[x]=(double)(ML[x]);
  85.      
  86.      if(M[0]==0.0)
  87.       {
  88.        ErrorHandler(OUT_OF_RANGE2,"Null mean in Classify");
  89.        M[0]=1E-100;
  90.       }
  91.      
  92.      mx=M[1]/M[0];
  93.      my=M[2]/M[0];
  94.      
  95.      CM[0] = M[0];              /* u_00 */
  96.      CM[1] = 0.0;               /* u_10 */
  97.      CM[2] = 0.0;               /* u_01 */
  98.      CM[3] = M[3]-my*M[1];      /* u_11 */
  99.      CM[4] = M[4]-mx*M[1];      /* u_20 */
  100.      CM[5] = M[5]-my*M[2];      /* u_02 */
  101.      CM[6] = M[6]-2.0*mx*M[3]-my*M[4]+2.0*mx*mx*M[2]; /* u_21 */
  102.      CM[7] = M[7]-2.0*my*M[3]-mx*M[5]+2.0*my*my*M[1]; /* u_12 */
  103.      CM[8] = M[8]-3.0*mx*M[4]+2.0*mx*mx*M[1];         /* u_30 */
  104.      CM[9] = M[9]-3.0*my*M[5]+2.0*my*my*M[2];         /* u_03 */
  105.      
  106.      iM2=1.0/(M[0]*M[0]);
  107.      iM3=1.0/pow(M[0],2.5);
  108.      
  109.      n[0]=CM[3]*iM2; /* n_11 */
  110.      n[1]=CM[4]*iM2; /* n_20 */
  111.      n[2]=CM[5]*iM2; /* n_02 */
  112.      n[3]=CM[6]*iM3; /* n_21 */
  113.      n[4]=CM[7]*iM3; /* n_12 */
  114.      n[5]=CM[8]*iM3; /* n_30 */
  115.      n[6]=CM[9]*iM3; /* n_03 */
  116.      
  117.      /* These features may give rise to numerical problems
  118.         with floats which for some systems must be in 
  119.         the range [1E-37;1E37]  */
  120.      
  121.      NCM[0]=n[1]+n[2];                               /* p_1 */
  122.      if (Dim>1) 
  123.       {
  124.        NCM[1]=Pow2(n[1]-n[2])+                       /* p_2 */
  125.        4.0*Pow2(n[0]);
  126.        if (Dim>2) 
  127.         {
  128.          NCM[2]=Pow2(n[5]-3.0*n[4])+                 /* p_3 */
  129.          Pow2(3.0*n[3]-n[6]);
  130.          if (Dim>3) 
  131.           {
  132.            NCM[3]=Pow2(n[5]+n[4])+                   /* p_4 */
  133.            Pow2(n[3]+n[6]);
  134.            if (Dim>4) 
  135.             {
  136.              NCM[4]=(n[5]-3.0*n[4])*(n[5]+n[4])*     /* p_5 */
  137.              (Pow2(n[5]+n[4])-
  138.               3.0*Pow2(n[3]+n[6]))+
  139.              (3.0*n[3]-n[6])*(n[3]+n[6])*
  140.              (3.0*Pow2(n[5]+n[4])-
  141.               Pow2(n[3]+n[6]));
  142.              if (Dim>5) 
  143.               {
  144.                NCM[5]=(n[1]-n[2])*                   /* p_6 */
  145.                (Pow2(n[5]+n[4])-      
  146.                 Pow2(n[3]+n[6]))+
  147.                4.0*n[0]*(n[5]+n[4])*
  148.                (n[3]+n[6]);
  149.                if (Dim>6) 
  150.                 {
  151.                  NCM[6]=(3.0*n[3]-n[5])*(n[5]+n[4])* /* p_7 */
  152.                  (Pow2(n[5]+n[4])-      
  153.                   3.0*Pow2(n[3]+n[6]))+
  154.                  (3.0*n[4]-n[5])*(n[3]+n[6])*
  155.                  (3.0*Pow2(n[5]+n[4])-
  156.                   Pow2(n[3]+n[6]));
  157.                 }
  158.               }
  159.             }
  160.           }
  161.         }
  162.       }
  163.      
  164.      for(x=0;x<Dim;x++)
  165.       {
  166.        temp=fabs(NCM[x]);
  167.        if (NCM[x]>=0.0) Features[x+PREDEFFEATURES]=(float)log(temp);
  168.        else             Features[x+PREDEFFEATURES]=(float)-log(temp);
  169.  
  170.        if (temp==0.0) 
  171.         {
  172.          vprintf(stderr,"Feature: %d",x); 
  173.          ErrorHandler(OUT_OF_RANGE2,"Feature value out of range");
  174.          Features[VAR]=-99999; /* replace with shade */
  175.         }
  176.       }
  177.      
  178.     }
  179.    else 
  180.     {
  181.      for(x=0;x<Dim;x++) Features[x+PREDEFFEATURES]=0.0;
  182.      Features[VAR] = -Features[VAR]; /* negative variance indicates a shade */
  183.     }
  184.    
  185.    return Features;
  186.   }
  187.